Effect of Cross-Linking Density on Non-Linear Viscoelasticity of Vulcanized SBR: A MD Simulation and Experimental Study

In recent years, there has been a growing interest in changes in dynamic mechanical properties of mixed rubber during dynamic shear, yet the influence of vulcanized characteristics on the dynamic shear behavior of vulcanized rubber, particularly the effect of cross-linking density, has received little attention. This study focuses on styrene–butadiene rubber (SBR) and aims to investigate the impact of different cross-linking densities (Dc) on dynamic shear behavior using molecular dynamics (MD) simulations. The results reveal a remarkable Payne effect, where the storage modulus experiences a significant drop when the strain amplitude (γ0) exceeds 0.1, which can be attributed to the fracture of the polymer bond and the decrease in the molecular chain’s flexibility. The influence of various Dc values mainly resides at the level of molecular aggregation in the system, where higher Dc values impede molecular chain motion and lead to an increase in the storage modulus of SBR. The MD simulation results are verified through comparisons with existing literature.


Introduction
Rubber is a cost-effective viscoelastic material with exceptional elasticity performance that finds extensive applications in automobiles, aerospace, construction, and medicine [1][2][3]. Since rubber products are usually subjected to alternating loads [4], the processing technology has a significant impact on their mechanical properties, particularly their dynamic mechanical properties. Studying the influence of processing technology on rubber's dynamic mechanical properties is crucial to enhancing its mechanical performance.
Pure rubber cannot meet practical mechanical performance requirements. Thus, the most common method to improve rubber mechanics is adding microfillers. Silicon dioxide particles and carbon black (CB) are the most common microfillers used to improve rubber mechanics. Jong et al. [5,6] showed that adding silica fillers significantly improves natural rubber's tensile properties, while Hait et al. [7] demonstrated the enhanced mechanical effect of adding CB to polybutadiene rubber (BR). Lin et al. [8] prepared CB-filled powdered natural rubber (P (NR/N234)) by the method of latex-CB coagulation technology and deeply studied the influence of curing recipes and CB contents on the curing properties, mechanical properties, and dynamic properties, and the results were compared with those of NR/N234 compounds based on traditional dry mixing of bale natural rubber and CB. Robertson and Hardman discussed in detail the mechanism of CB reinforcing the mechanical properties of rubber and the generation mechanism of the Payne effect in the shear process of filling rubbers [9], and there is no doubt that this work is helpful for understanding the reinforcing mechanism and Payne effect of fillers on rubber.
The Payne effect relates to the elastic and storage moduli variation, which typically increase and decrease, respectively, under high shear strain. The phenomenon was first reported by Payne and has since garnered considerable attention. It has been extensively applied in dynamic mechanical analysis for polymer materials under medium-to-high strain. The addition of nanoparticles (NPs) to the mixed rubber has a significant impact on the Payne effect. Shi et al.'s research studied the enhancement mechanism of the Payne effect by using CB to fill natural rubber, and the results indicated the recoverability and hysteresis of the Payne effect, demonstrating the positive impact of high frequencies and low temperatures on enhanced Payne effect. The structural evolution of the filler phase was not the primary influence factor [10]. Similarly, Zhao et al. [11] applied CB as a filler to SBR to create a CB conductive network and studied the Payne effect of the filling system. The results revealed that the polymer matrix next to CB particles can accelerate the reconstruction of the fillerpolymer network, thereby enhancing the reversibility of the Payne effect [12][13][14]. CB, along with other fillers, such as spherical [15][16][17] and fibrous [18][19][20], have excellent electrical and thermal conductivity, which significantly affects Payne effect improvement. Experimental research can satisfactorily explain the Payne effect of mixed rubber filling. However, the micromechanism of Payne effect formation remains challenging to explain.
Molecular dynamics (MD) simulation is often preferred to experiments for studying the microscopic dynamic mechanisms of the Payne effect in rubber. MD has several advantages over experiments because it can explore the impact of molecular chain motion and filler-polymer interaction forces, which are typically unknown under experimental conditions. Hong et al. [21] demonstrated the utility of non-equilibrium molecular dynamics (NEMD) in exploring the viscoelastic properties of nanoplate-filled polymer composites as it relates to filler packing fraction, filler-polymer interaction forces, Rouse dynamics, chain constraint, and percolation networks. Gao et al. [22] used hydrodynamic effects, "clog rubber" and "filler network", to explain how polymer nanocomposites filled with nanorods experience dynamic modulus reinforcement and showed that the Payne effect becomes more pronounced with stronger interface interaction and increased nanorod fraction.
Currently, most research on the Payne effect has focused on mixed rubber, with little research done on vulcanized rubber. Vulcanization plays a critical role in rubber processing, but existing research about the impact of vulcanization characteristics on the Payne effect of rubber is limited. Therefore, it is necessary to investigate the effect of different vulcanization characteristics on the dynamic mechanical properties of rubber. Hou et al. [23] conducted experiments on the Payne effect of polyisoprene rubber vulcanized using thermooxidative aging, and the results revealed that thermooxidative aging significantly affected the vulcanized structure and Payne effect of polyisoprene rubber. Employing a highly efficient vulcanization method to prepare vulcanized polyisoprene rubber will remarkably enhance Payne dissipation, as for its microscopic mechanism, it was not mentioned in this report.
The application scope of vulcanized SBR is wide; therefore, the study of its Payne effect holds theoretical and practical significance. This research aims to utilize molecular dynamics methods to explore the Payne effect of vulcanized SBR in dynamic shear processes, providing new ideas for the study of the dynamic mechanical behavior of sulfurized rubber materials.

Verification of Simulation Correctness
The cross-linking densities (mol/mL) corresponding to each D c in MD model obtained via calculation are presented in Table 1. The results of cross-linking densities obtained by experiments are also shown in Table 1.
The small difference between calculations and the experimental results suggests that MD simulation is able to realistically simulate the experiment, thereby demonstrating high reliability in this study.
In order to prove the accuracy of simulation model, we perform a comparison between T g of SBR systems exhibiting varying D c , utilizing both molecular simulation and experimentation approaches. Figure 1a-e depict the relationships between density and temperature of SBR systems at different D c through molecular dynamic simulation. Figure 1f illustrates the T g of each sample by experiments. Table 2 enumerates the T g for SBR systems of different D c acquired through both the experimental and molecular dynamic simulation methods, alongside the corresponding relative deviation. By comparing the DSC results of different S contents with the T g results obtained through simulation (Figure 1), it can be found that the results of MD simulation can predict the T g relatively well. In order to prove the accuracy of simulation model, we perform a comparison between Tg of SBR systems exhibiting varying Dc, utilizing both molecular simulation and experimentation approaches. Figure 1a-e depict the relationships between density and temperature of SBR systems at different Dc through molecular dynamic simulation. Figure  1f illustrates the Tg of each sample by experiments. Table 2 enumerates the Tg for SBR systems of different Dc acquired through both the experimental and molecular dynamic simulation methods, alongside the corresponding relative deviation. By comparing the DSC results of different S contents with the Tg results obtained through simulation (Figure 1), it can be found that the results of MD simulation can predict the Tg relatively well.  The tensile curve of pure SBR at strain rate of 0.01 ps −1 is shown in Figure 2, where snapshots captured at strains of 0, 1, and 2 are also provided. Additionally, Figure 2 presents the tensile curve obtained using the identical simulation parameters in reference [24]. It can be found in Figure 2 that the tensile strength of pure SBR is approximately 2 GPa at a strain rate of 0.01 ps −1 , consistent with the data reported in [24]. The high reliability of the pcff force field we selected is thus demonstrated.

The Effects of Dc and γ0 on Dynamic Mechanical Properties
The stress-strain curve of SBR with Dc of 8.0 at γ0 = 0.5 is extracted to observe the phenomenon of strain lagging behind stress. Figure 3 presents the stress-strain curve, which shows the strain curve is constantly delayed behind the stress curve by a constant time δt of 20 fs. Notably, the stress curve's first two peaks are higher than the others. The phenomenon is predominantly attributed to the hindered movement of the entangled polymer chains by an external force, i.e., the shear stress. The interference results in relatively high viscosity and stress, leading to slightly higher peak values for the first two peaks. The originally entangled polymer chains' entanglement points are released with the continued shearing, reducing the viscosity and leading to the gradual stabilization of the stress curve's peak values. The tensile curve of pure SBR at strain rate of 0.01 ps −1 is shown in Figure 2, where snapshots captured at strains of 0, 1, and 2 are also provided. Additionally, Figure 2 presents the tensile curve obtained using the identical simulation parameters in reference [24]. It can be found in Figure 2 that the tensile strength of pure SBR is approximately 2 GPa at a strain rate of 0.01 ps −1 , consistent with the data reported in [24]. The high reliability of the pcff force field we selected is thus demonstrated.
. The tensile curve of pure SBR at strain rate of 0.01 ps −1 is shown in Figure 2, where snapshots captured at strains of 0, 1, and 2 are also provided. Additionally, Figure 2 presents the tensile curve obtained using the identical simulation parameters in reference [24]. It can be found in Figure 2 that the tensile strength of pure SBR is approximately 2 GPa at a strain rate of 0.01 ps −1 , consistent with the data reported in [24]. The high reliability of the pcff force field we selected is thus demonstrated.

The Effects of Dc and γ0 on Dynamic Mechanical Properties
The stress-strain curve of SBR with Dc of 8.0 at γ0 = 0.5 is extracted to observe the phenomenon of strain lagging behind stress. Figure 3 presents the stress-strain curve, which shows the strain curve is constantly delayed behind the stress curve by a constant time δt of 20 fs. Notably, the stress curve's first two peaks are higher than the others. The phenomenon is predominantly attributed to the hindered movement of the entangled polymer chains by an external force, i.e., the shear stress. The interference results in relatively high viscosity and stress, leading to slightly higher peak values for the first two peaks. The originally entangled polymer chains' entanglement points are released with the continued shearing, reducing the viscosity and leading to the gradual stabilization of the stress curve's peak values.

The Effects of D c and γ 0 on Dynamic Mechanical Properties
The stress-strain curve of SBR with D c of 8.0 at γ 0 = 0.5 is extracted to observe the phenomenon of strain lagging behind stress. Figure 3 presents the stress-strain curve, which shows the strain curve is constantly delayed behind the stress curve by a constant time δ t of 20 fs. Notably, the stress curve's first two peaks are higher than the others. The phenomenon is predominantly attributed to the hindered movement of the entangled polymer chains by an external force, i.e., the shear stress. The interference results in relatively high viscosity and stress, leading to slightly higher peak values for the first two peaks. The originally entangled polymer chains' entanglement points are released with the continued shearing, reducing the viscosity and leading to the gradual stabilization of the stress curve's peak values.   Figure 4 illustrates that SBR with varying Dc manifests conspicuous Payne effects under varying γ0. Specifically, the G′ experiences a drastic decline with the augmentation of γ0, as depicted in Figure 4a. Additionally, Figure 4a shows that the G′ escalates with the intensification of Dc. Figure 4b indicates that the tanδ of different Dc displays an increasing tendency as γ0 increases. To explore the formation mechanism of Payne effects and the Dc impact on dynamic performance during the shear process, detailed analysis from both a static and dynamic outlook are operated.

Static Equilibrium Process
MSD and RDF of C atoms on the main chain are extracted from Dc systems with differing properties at equilibrium at the NPT for a period of 2 ns and presented in Figures  5 and 6. The purpose of this extraction is to investigate the impact of movement of molecular chains and arrangement of atomic components on mechanical characteristics.    Figure 4 illustrates that SBR with varying Dc manifests conspicuous Payne effects under varying γ0. Specifically, the G′ experiences a drastic decline with the augmentation of γ0, as depicted in Figure 4a. Additionally, Figure 4a shows that the G′ escalates with the intensification of Dc. Figure 4b indicates that the tanδ of different Dc displays an increasing tendency as γ0 increases. To explore the formation mechanism of Payne effects and the Dc impact on dynamic performance during the shear process, detailed analysis from both a static and dynamic outlook are operated.

Static Equilibrium Process
MSD and RDF of C atoms on the main chain are extracted from Dc systems with differing properties at equilibrium at the NPT for a period of 2 ns and presented in Figures  5 and 6. The purpose of this extraction is to investigate the impact of movement of molecular chains and arrangement of atomic components on mechanical characteristics.  Figure 4 illustrates that SBR with varying D c manifests conspicuous Payne effects under varying γ 0 . Specifically, the G experiences a drastic decline with the augmentation of γ 0 , as depicted in Figure 4a. Additionally, Figure 4a shows that the G escalates with the intensification of D c . Figure 4b indicates that the tanδ of different D c displays an increasing tendency as γ 0 increases. To explore the formation mechanism of Payne effects and the D c impact on dynamic performance during the shear process, detailed analysis from both a static and dynamic outlook are operated.

Static Equilibrium Process
MSD and RDF of C atoms on the main chain are extracted from D c systems with differing properties at equilibrium at the NPT for a period of 2 ns and presented in Figures 5 and 6. The purpose of this extraction is to investigate the impact of movement of molecular chains and arrangement of atomic components on mechanical characteristics. Observations from Figure 5 reveal that the MSD exhibits a decreasing trend with an increasing Dc. Upon attaining equilibrium, the MSD of Dc = 1.0 SBR system is roughly 6 times more considerable than that of the Dc = 8.0 system. This is due to the abundance of C-S-C bonds, which progressively limit the movement of molecular chains, restricting their motion range and consequently their rigidity. This phenomenon, in turn, explains the increasing trend in G′ with Dc increment, as illustrated in Figure 4a. The RDF depicted in Figure 6 reveals an increased density of C atoms present on the SBR main chain with an increment of Dc. Specifically, for the Dc = 8.0 system, this density level is observed to be 10% higher than for the Dc = 1.0 counterpart. The reason for this occurrence is primarily attributed to the presence of more C-S-C bonds that support a tighter coupling between the SBR chains, as depicted in Figure 5.

Microscopic Mechanism of Payne Effect
In this work, we sought to examine the influence of cross-links on the shearing process by calculating the rate of broken cross-linking bonds, or P, after10 cycles of shearing. The calculation is performed as follows. If the length of a C-S bond linked with an S atom exceeded 1.8 Å, it is designated as "broken" and considered unrepairable during the shearing process. P is calculated from the ratio of broken cross-links (Nb) to the total number of cross-linking bonds (Ns), as given by P = Nb/Ns, where the subscript "b" represents Observations from Figure 5 reveal that the MSD exhibits a decreasing trend with an increasing D c . Upon attaining equilibrium, the MSD of D c = 1.0 SBR system is roughly 6 times more considerable than that of the D c = 8.0 system. This is due to the abundance of C-S-C bonds, which progressively limit the movement of molecular chains, restricting their motion range and consequently their rigidity. This phenomenon, in turn, explains the increasing trend in G with D c increment, as illustrated in Figure 4a. Observations from Figure 5 reveal that the MSD exhibits a decreasing trend with an increasing Dc. Upon attaining equilibrium, the MSD of Dc = 1.0 SBR system is roughly 6 times more considerable than that of the Dc = 8.0 system. This is due to the abundance of C-S-C bonds, which progressively limit the movement of molecular chains, restricting their motion range and consequently their rigidity. This phenomenon, in turn, explains the increasing trend in G′ with Dc increment, as illustrated in Figure 4a. The RDF depicted in Figure 6 reveals an increased density of C atoms present on the SBR main chain with an increment of Dc. Specifically, for the Dc = 8.0 system, this density level is observed to be 10% higher than for the Dc = 1.0 counterpart. The reason for this occurrence is primarily attributed to the presence of more C-S-C bonds that support a tighter coupling between the SBR chains, as depicted in Figure 5.

Microscopic Mechanism of Payne Effect
In this work, we sought to examine the influence of cross-links on the shearing process by calculating the rate of broken cross-linking bonds, or P, after10 cycles of shearing. The calculation is performed as follows. If the length of a C-S bond linked with an S atom exceeded 1.8 Å, it is designated as "broken" and considered unrepairable during the shearing process. P is calculated from the ratio of broken cross-links (Nb) to the total number of cross-linking bonds (Ns), as given by P = Nb/Ns, where the subscript "b" represents The RDF depicted in Figure 6 reveals an increased density of C atoms present on the SBR main chain with an increment of D c . Specifically, for the D c = 8.0 system, this density level is observed to be 10% higher than for the D c = 1.0 counterpart. The reason for this occurrence is primarily attributed to the presence of more C-S-C bonds that support a tighter coupling between the SBR chains, as depicted in Figure 5.

Microscopic Mechanism of Payne Effect
In this work, we sought to examine the influence of cross-links on the shearing process by calculating the rate of broken cross-linking bonds, or P, after10 cycles of shearing. The calculation is performed as follows. If the length of a C-S bond linked with an S atom exceeded 1.8 Å, it is designated as "broken" and considered unrepairable during the shearing process. P is calculated from the ratio of broken cross-links (N b ) to the total number of cross-linking bonds (N s ), as given by P = N b /N s , where the subscript "b" represents "broken". We present the relationship between P and γ 0 for various SBR systems with differing D c values in Figure 7.
"broken". We present the relationship between P and γ0 for various SBR systems with differing Dc values in Figure 7. By analyzing Figure 7, it is evident that an increase in γ0 value leads to a corresponding increase in P. Larger γ0 values in a vulcanized SBR system cause the breakage of C-S-C bonds, which increases as γ0 increases. When γ0 exceeds 0.1, all SBR systems demonstrate a sharp rise in P value, consistent with the observation in Figure 4a. Therefore, it can be concluded that larger values of γ0 lead to a sudden increase in the number of broken cross-linking bonds, resulting in the occurrence of Payne effects during the shearing process, affecting the mechanical properties of the SBR system. Furthermore, the P value decreases as Dc increases. This implies that higher cross-linked densities exhibit lower P values after shearing. A plausible explanation for this observation is that larger Dc values can enhance the overall mechanical performance of SBR. During the shearing process, systems with lower Dc values tend to break more easily than those with higher Dc values, resulting in lower P values. The breakage of cross-linking bonds during the shearing process is mainly related to the movement of molecular chains. Therefore, taking the system with Dc = 8.0 as an example, the MSD during the shearing process is extracted under different values of γ0, as shown in Figure 8. In addition, Figure 8 also includes the MSD of different SBR systems with γ0 = 0.5 and varying Dc. By analyzing Figure 7, it is evident that an increase in γ 0 value leads to a corresponding increase in P. Larger γ 0 values in a vulcanized SBR system cause the breakage of C-S-C bonds, which increases as γ 0 increases. When γ 0 exceeds 0.1, all SBR systems demonstrate a sharp rise in P value, consistent with the observation in Figure 4a. Therefore, it can be concluded that larger values of γ 0 lead to a sudden increase in the number of broken crosslinking bonds, resulting in the occurrence of Payne effects during the shearing process, affecting the mechanical properties of the SBR system. Furthermore, the P value decreases as D c increases. This implies that higher cross-linked densities exhibit lower P values after shearing. A plausible explanation for this observation is that larger D c values can enhance the overall mechanical performance of SBR. During the shearing process, systems with lower D c values tend to break more easily than those with higher D c values, resulting in lower P values. The breakage of cross-linking bonds during the shearing process is mainly related to the movement of molecular chains. Therefore, taking the system with D c = 8.0 as an example, the MSD during the shearing process is extracted under different values of γ 0 , as shown in Figure 8. In addition, Figure 8 also includes the MSD of different SBR systems with γ 0 = 0.5 and varying D c . "broken". We present the relationship between P and γ0 for various SBR systems with differing Dc values in Figure 7. By analyzing Figure 7, it is evident that an increase in γ0 value leads to a corresponding increase in P. Larger γ0 values in a vulcanized SBR system cause the breakage of C-S-C bonds, which increases as γ0 increases. When γ0 exceeds 0.1, all SBR systems demonstrate a sharp rise in P value, consistent with the observation in Figure 4a. Therefore, it can be concluded that larger values of γ0 lead to a sudden increase in the number of broken cross-linking bonds, resulting in the occurrence of Payne effects during the shearing process, affecting the mechanical properties of the SBR system. Furthermore, the P value decreases as Dc increases. This implies that higher cross-linked densities exhibit lower P values after shearing. A plausible explanation for this observation is that larger Dc values can enhance the overall mechanical performance of SBR. During the shearing process, systems with lower Dc values tend to break more easily than those with higher Dc values, resulting in lower P values. The breakage of cross-linking bonds during the shearing process is mainly related to the movement of molecular chains. Therefore, taking the system with Dc = 8.0 as an example, the MSD during the shearing process is extracted under different values of γ0, as shown in Figure 8. In addition, Figure 8 also includes the MSD of different SBR systems with γ0 = 0.5 and varying Dc.  It can be observed in Figure 8a that MSD is approximately 7.0 at γ0 = 0.09, but it significantly increases when γ0 exceeds 0.1, as seen in Figure 8e-g. For example, when γ0 = 0.5, MSD reaches approximately 230, which is 33 times that of γ0 = 0.09. Higher γ0 values move the system away from equilibrium, increasing the probability of cross-linking bond breakage and potentially decreasing mechanical performance, as explained by the significant decline in G′ shown in Figure 4a when γ0 is greater than 0.1. Furthermore, MSD decreases as Dc increases under the same γ0 value, as shown in Figure 8c,d. This is primarily due to the partial restriction of molecular chain movement by the formation of the crosslinking network, consequently slowing the occurrence of bond breakage. Hence, when γ0 remains the same, an increase in Dc results in a higher G′.
To analyze the energy variation in SBR systems during the shearing process, we calculate both the non-bonding energy (Epair, which describes the repulsive or attractive forces between two non-bonded atoms and can be calculated by the sum of the van der Waals energy and Coulomb energy) as well as the bonding energy, which included bond energy (Ebond, which represents the stretching energy of a bond and characterizes the energy change resulting from the movement of each chemical bond in the direction of the molecular axis), angle energy (Eangle, which represents the bending energy of a bond angle, caused by changes to the bond angle) and dihedral energy (Edihed, which represents the energy change caused by the rotation of a single bond and resulting in molecular skeleton distortion) over the period of shearing. Initially, we extracted Etotal, Epair, Ebond, Eangle, and Edihed of SBR systems with different values of Dc during the shearing process at γ0 = 0.5, as depicted in Figure 9 (where Etotal = Epair + Ebond + Eangle + Edihed). We also present the energy variation of a non-cross-linking system SBR system (Dc = 0) during the shearing process in Figure 9 for comparing the specific effects of different cross-linking densities on energy changes. We consider the energy difference (E-E0) at each moment of the shearing process relative to the system's initial energy [25] to facilitate analysis. It can be observed in Figure 8a that MSD is approximately 7.0 at γ 0 = 0.09, but it significantly increases when γ 0 exceeds 0.1, as seen in Figure 8e-g. For example, when γ 0 = 0.5, MSD reaches approximately 230, which is 33 times that of γ 0 = 0.09. Higher γ 0 values move the system away from equilibrium, increasing the probability of cross-linking bond breakage and potentially decreasing mechanical performance, as explained by the significant decline in G shown in Figure 4a when γ 0 is greater than 0.1. Furthermore, MSD decreases as D c increases under the same γ 0 value, as shown in Figure 8c,d. This is primarily due to the partial restriction of molecular chain movement by the formation of the cross-linking network, consequently slowing the occurrence of bond breakage. Hence, when γ 0 remains the same, an increase in D c results in a higher G .
To analyze the energy variation in SBR systems during the shearing process, we calculate both the non-bonding energy (E pair , which describes the repulsive or attractive forces between two non-bonded atoms and can be calculated by the sum of the van der Waals energy and Coulomb energy) as well as the bonding energy, which included bond energy (E bond , which represents the stretching energy of a bond and characterizes the energy change resulting from the movement of each chemical bond in the direction of the molecular axis), angle energy (E angle , which represents the bending energy of a bond angle, caused by changes to the bond angle) and dihedral energy (E dihed , which represents the energy change caused by the rotation of a single bond and resulting in molecular skeleton distortion) over the period of shearing. Initially, we extracted E total , E pair , E bond , E angle , and E dihed of SBR systems with different values of D c during the shearing process at γ 0 = 0.5, as depicted in Figure 9 (where E total = E pair + E bond + E angle + E dihed ). We also present the energy variation of a non-cross-linking system SBR system (D c = 0) during the shearing process in Figure 9 for comparing the specific effects of different cross-linking densities on energy changes. We consider the energy difference (E-E 0 ) at each moment of the shearing process relative to the system's initial energy [25] to facilitate analysis. It is obvious in Figure 9 that the energy of the various Dc systems fluctuates si dally during shear and gradually decreases. This behavior can be adequately exp by combining Figure 3b which results from the peak stress increase. Additionally, crease in Dc leads to a rise in Etotal, which increases by approximately 12% for D compared to the non-cross-linking system. This can mainly be attributed to the fa the higher Dc system has more cross-linking bonds, thus increasing the Etotal of the sy Furthermore, it can be observed from Figure 9b that Epair does not significantly cha Dc increases. This is due to SBR chains being mainly cross-linked by C-S-C bonds, me that interatomic interaction forces do not show significant differences.
In contrast to Epair, both Ebond and Eangle corresponding to Figure 9c,d exhibit corre with Dc. The peak values of Ebond and Eangle increase and show an increase of 19.4% and respectively as Dc increases, when compared to the non-cross-linking system. Fig  shows that there is no significant correlation with respect to changes in Dc for Edihed. It is obvious in Figure 9 that the energy of the various D c systems fluctuates sinusoidally during shear and gradually decreases. This behavior can be adequately explained by combining Figure 3b which results from the peak stress increase. Additionally, an increase in D c leads to a rise in E total , which increases by approximately 12% for D c = 8.0 compared to the non-cross-linking system. This can mainly be attributed to the fact that the higher D c system has more cross-linking bonds, thus increasing the E total of the system. Furthermore, it can be observed from Figure 9b that E pair does not significantly change as D c increases. This is due to SBR chains being mainly cross-linked by C-S-C bonds, meaning that interatomic interaction forces do not show significant differences.
In contrast to E pair , both E bond and E angle corresponding to Figure 9c,d exhibit correlation with D c . The peak values of E bond and E angle increase and show an increase of 19.4% and 6.5%, respectively as D c increases, when compared to the non-cross-linking system. Figure 9e shows that there is no significant correlation with respect to changes in D c for E dihed .
To sum up, from an energy perspective, changing D c mainly affects the stretching of bonds and the rotation of bond angles during the shearing process under the same γ 0 conditions. As for the rotation of bond angles, D c has a more significant impact on bond stretching in the cross-linking SBR system. In order to investigate the impact of different γ 0 on the energy during the shearing process, we extract the energy variation of the system with D c = 8.0 under different γ 0 , as illustrated in Figure 10.
To sum up, from an energy perspective, changing Dc mainly affects the stretching of bonds and the rotation of bond angles during the shearing process under the same γ0 conditions. As for the rotation of bond angles, Dc has a more significant impact on bond stretching in the cross-linking SBR system. In order to investigate the impact of different γ0 on the energy during the shearing process, we extract the energy variation of the system with Dc = 8.0 under different γ0, as illustrated in Figure 10. According to Figure 10a, the Etotal displays periodic variation with two peaks during each shear cycle, corresponding to the maximal and minimal strains of the action, respectively. As the parameter γ0 increases, the peaks of Etotal also increase. This finding has been previously discussed. The process of shearing depletes Etotal gradually, primarily due to According to Figure 10a, the E total displays periodic variation with two peaks during each shear cycle, corresponding to the maximal and minimal strains of the action, respectively. As the parameter γ 0 increases, the peaks of E total also increase. This finding has been previously discussed. The process of shearing depletes E total gradually, primarily due to the viscoelasticity of SBR absorbing energy through overcoming the frictional force of molecu-lar chains [26]. Consequently, the energy reduction rate tends to be slower under lower strain conditions, when molecular movement gets more limited, and resulting in lower frictional resistance between chains, causing slower dissipation of energy. Figure 10c,d report alterations in E bond and E angle during shearing. The shift trends of these values lead to the conclusion that the energy change during shearing is mainly driven by bond stretching and angle variation. The sum of E bond and E angle accounts for more than half of E total , making these the primary contributors to the shearing status of the system. During the shearing process, stretching bonds and rotating bond angles mainly determine the deformation mode of the system. E bond tends to reach equilibrium at 20,000 kcal/mol and E angle at 14,000 kcal/mol. These findings suggest that bond stretching has a superior role concerning bond angles rotation in the shearing process. Dihedral angle distortion plays a relatively minor part in contrast, Figure 10e shows that E dihed exhibits periodic variation and increases with γ 0 , while gradually decreasing with an increase in shearing cycle. When γ 0 = 0.05, E dihed is around −400 kcal/mol and when γ 0 = 0.5, E dihed is around 4000 kcal/mol. The principal reason for the raised value of E dihed is the increase in the distance between the molecular chains and their initial position due to a larger γ 0 .
To study the changes in molecular chain conformation during the shearing process, we extract the mean square radius of gyration (R g ) of the molecular chain at D c = 8.0 under different γ 0 , calculated by Equation (1) and shown in Figure 11.
where M is the total of the model, r cm is the center-of-mass position, and the sum is over all atoms in model. The T g is the primary mechanism for examining the impact of molecular chain flexibility on dynamic mechanical properties. A larger R g indicates poorer flexibility and molecular chains with poorer flexibility correspond to higher T g and poorer dynamic mechanical performance at the same temperature [27,28].
From Figure 11a-c, it is evident that the R g of the molecular chain undergoes periodic variation during the shearing process, with two peaks in one shear cycle, corresponding to the two displacement limits during the shearing process. With the increase in parameter γ 0 , the R g peaks increase, implying a decreased flexibility of molecular chains during shearing, ultimately leading to diminished dynamic mechanical properties of vulcanized SBR. the viscoelasticity of SBR absorbing energy through overcoming the frictional force of molecular chains [26]. Consequently, the energy reduction rate tends to be slower under lower strain conditions, when molecular movement gets more limited, and resulting in lower frictional resistance between chains, causing slower dissipation of energy. Figure  10c,d report alterations in Ebond and Eangle during shearing. The shift trends of these values lead to the conclusion that the energy change during shearing is mainly driven by bond stretching and angle variation. The sum of Ebond and Eangle accounts for more than half of Etotal, making these the primary contributors to the shearing status of the system. During the shearing process, stretching bonds and rotating bond angles mainly determine the deformation mode of the system. Ebond tends to reach equilibrium at 20,000 kcal/mol and Eangle at 14,000 kcal/mol. These findings suggest that bond stretching has a superior role concerning bond angles rotation in the shearing process. Dihedral angle distortion plays a relatively minor part in contrast, Figure 10e shows that Edihed exhibits periodic variation and increases with γ0, while gradually decreasing with an increase in shearing cycle. When γ0 = 0.05, Edihed is around −400 kcal/mol and when γ0 = 0.5, Edihed is around 4000 kcal/mol. The principal reason for the raised value of Edihed is the increase in the distance between the molecular chains and their initial position due to a larger γ0.
To study the changes in molecular chain conformation during the shearing process, we extract the mean square radius of gyration (Rg) of the molecular chain at Dc = 8.0 under different γ0, calculated by Equation (1) and shown in Figure 11.
where M is the total of the model, rcm is the center-of-mass position, and the sum is over all atoms in model. The Tg is the primary mechanism for examining the impact of molecular chain flexibility on dynamic mechanical properties. A larger Rg indicates poorer flexibility and molecular chains with poorer flexibility correspond to higher Tg and poorer dynamic mechanical performance at the same temperature [27,28].  From Figure 11a-c, it is evident that the Rg of the molecular chain undergoes periodic variation during the shearing process, with two peaks in one shear cycle, corresponding to the two displacement limits during the shearing process. With the increase in parameter γ0, the Rg peaks increase, implying a decreased flexibility of molecular chains during shearing, ultimately leading to diminished dynamic mechanical properties of vulcanized SBR.
The results from Figure 11 reveal that Rg remains almost unchanged, i.e., independent of the γ0 when γ0 is less than 0.1 (Figure 11c). On the contrary, Rg demonstrates a noticeable upward pattern with γ0, as determined by the difference between the initial Rg and the value after 10 shearing cycles, when γ0 exceeds 0.1 (Figure 11b,d). As the value of γ0 increases, the molecular chains become progressively inflexible, causing a significant increase in Rg and a reduction in the mechanical properties. This trend is further demonstrated in Figure 4a, which illustrates the G′ of SBR, showing no obvious reduction when γ0 falls below 0.1. However, when γ0 reaches and surpasses 0.1, the G′ of SBR dwindles considerably.

Verification of Simulation Correctness
The cross-linking densities (mol/mL) corresponding to each Dc in MD model obtained via calculation are presented in Table 3. The results of cross-linking densities obtained by experiments are also shown in Table 3.
where Ns is the number of cross-link bonds and Mc is the number of single chains. The cross-linking criteria are defined as follows. If the distance between two carbon atoms (C) on the main chain falls within the range of 1.6 to 3.45 Å, a sulfur atom (S) is added to the geometric center of the two C atoms and a C-S-C bond is formed [29]. Additionally, since every two C atoms on a single chain can create a C-S-C bond, the phenomenon of self-cross-linking is considered in the actual cross-linking process [30]. Once a C atom on the main chain establishes a cross-link bond, it can no longer bind to other S Figure 11. R g of the system under different γ 0 (a-c) varies with the shear period and the difference between the R g and the initial R g ((R g ) T -(R g ) 0 ) after 10 shear periods under different γ 0 (d).
The results from Figure 11 reveal that R g remains almost unchanged, i.e., independent of the γ 0 when γ 0 is less than 0.1 (Figure 11c). On the contrary, R g demonstrates a noticeable upward pattern with γ 0 , as determined by the difference between the initial R g and the value after 10 shearing cycles, when γ 0 exceeds 0.1 (Figure 11b,d). As the value of γ 0 increases, the molecular chains become progressively inflexible, causing a significant increase in R g and a reduction in the mechanical properties. This trend is further demonstrated in Figure 4a, which illustrates the G of SBR, showing no obvious reduction when γ 0 falls below 0.1. However, when γ 0 reaches and surpasses 0.1, the G of SBR dwindles considerably.

Verification of Simulation Correctness
The cross-linking densities (mol/mL) corresponding to each D c in MD model obtained via calculation are presented in Table 3. The results of cross-linking densities obtained by experiments are also shown in Table 3. Table 3. Monomer components in single-chain SBR.

Monomer
where Ns is the number of cross-link bonds and Mc is the number of single chains. The cross-linking criteria are defined as follows. If the distance between two carbon atoms (C) on the main chain falls within the range of 1.6 to 3.45 Å, a sulfur atom (S) is added to the geometric center of the two C atoms and a C-S-C bond is formed [29]. Additionally, since every two C atoms on a single chain can create a C-S-C bond, the phenomenon of self-cross-linking is considered in the actual cross-linking process [30]. Once a C atom on the main chain establishes a cross-link bond, it can no longer bind to other S atoms. To implement the complete cross-linking process, this study utilized Perl scripts inside Materials Studio software.
The pcff force field is utilized to describe the intermolecular and intramolecular interactions in SBR. Equation (3) depicts the formula of the pcff force field.
where ε is the bond energy at the equilibrium position where the force F(r) = 0, σ is the collision diameter, r c is the truncation radius of 9.5 Å, and r refers to the interatomic distance. Figure 12 presents the molecular formula, process and results of cross-linking. After the completion of cross-linking, we perform relaxation annealing on the model. The specific method is as follows. Firstly, a 2 ns NPT relaxation is carried out at 300 K and 1 atm, then the model is heated to 600 K and then cooled to 300 K in NVT ensembles to simulate the annealing process. This process lasted for 4 ns and is repeated 4 times. Finally, a 2 ns NPT relaxation is performed at 300 K and 1 atm, during which the radial distribution function (RDF) and mean square displacement (MSD) of the system are extracted. The obtained model is used for studying the dynamic shear process.
This study investigated the dynamic mechanical properties of SBR with varying density of cross-linking (D c ) through cyclic shear deformation of simulated box. The shear deformation is performed in the XY plane by moving the simulation box along the ±X directions with a fixed shear rate of 0.01 ps −1 . The shear strain amplitude (γ 0 ) varied from 0.05 to 0.5 and the system underwent 10 complete shear cycles. The changes over time in shear stress (σ) and strain (γ) are expressed using Equations (5) and (6). Consequently, the storage modulus (G ) and loss factor (tanδ) of SBR with varying D c are obtained through Equation (7).
An Anderson thermostat and Berendsen barostat are employed to control temperature and pressure, respectively, for all simulation in this study, while the velocity Verlet algorithm is used for kinetic integration. Molecular dynamics simulations are performed This study investigated the dynamic mechanical properties of SBR with varying density of cross-linking (Dc) through cyclic shear deformation of simulated box. The shear deformation is performed in the XY plane by moving the simulation box along the ±X directions with a fixed shear rate of 0.01 ps −1 . The shear strain amplitude (γ0) varied from 0.05 to 0.5 and the system underwent 10 complete shear cycles. The changes over time in shear stress (σ) and strain (γ) are expressed using Equations (5) and (6). Consequently, the storage modulus (G′) and loss factor (tanδ) of SBR with varying Dc are obtained through Equation (7).
An Anderson thermostat and Berendsen barostat are employed to control temperature and pressure, respectively, for all simulation in this study, while the velocity Verlet algorithm is used for kinetic integration. Molecular dynamics simulations are performed using the LAMMPS package on a 24-core supercomputer [31]. The visualizing results are analyzed through Open Visualization Tool (OVITO) software [32].

Experimental Methods
The amounts of sulfur addition in vulcanized SBR are determined based on the S atomic mass fraction (phr) in different D c models. The determined amounts are 1.04, 1.99, 5.03, 6.29, and 8.28, respectively. The mixed SBR is obtained by adding compounding agents to pure SBR. Table 4 illustrates the content of each component. Pure SBR is premixed for 2 min on an open mill with a diameter of Φ152.4 mm. Subsequently, ZnO, SA, accelerant DZ, accelerant TT, and insoluble sulfur (S) are added in sequence to produce mixed SBR. The mixed rubber is allowed to stand for 24 h before testing the vulcanization properties and recording T 90 using a non-rotor vulcanizer produced by Gotech. Vulcanizations of each mixed rubber are carried out at 150 • C to create a 2 mm vulcanized SBR sheet with a temperature of 150 • C, and the time of vulcanization is calculated based on T 90 value.

Cross-Linking Density Test
After allowing the vulcanized SBR to stand for 24 h, the cross-linking density analyzer with the model XLDS-15 is used to test the cross-linking density by the method of nuclear magnetic resonance (NMR). The SBR samples are tested under specific conditions of 90 • C for the temperature, 0.35 T for the magnetic field intensity and 15 MHz for the resonance frequency. The cross-linking density is calculated by measuring the relaxation time (T 2 ) of the SBR samples with varying sulfur contents. The same formulae are tested five times and the average values are adopted as the final result.

DSC Test
The differential scanning calorimeter model DSC 200 F3 is utilized to quantify the glass transition temperature (T g ) of the vulcanized SBR. The sample weighing 8-10 mg is then loaded into the calorimeter and heated to 100 • C for 5 min to eliminate the thermal history. Next, the temperature is lowered to −75 • C and held for five minutes before being heated back to 100 • C. The test parameters include both the heating and cooling rates of 10 • C/min.

RPA Test
The RPA2000 rubber processing analyzer (RPA) is employed for evaluation of the dynamic mechanical properties of each sample. Initially, the temperature is set to 150 • C and the samples are subsequently put into the testing equipment while recording the time taken. Upon reaching the corresponding T 90 to simulate the vulcanization process, the instruments' temperature is decreased to 20 • C, after which the samples undergo a shearing test, the test frequency is set at 10 Hz and strain range is 0-55%.  Figure 13 shows an increase in the Payne effect in all groups of SBR with an increase in γ0. As observed from combining these findings with the results presented in Figure 5, it is apparent that higher S content enhances the deformation resistance of SBR, which is demonstrated by a noticeable upward trend in G′ ( Figure 13). Furthermore, as γ0 exceeds 0.1, G′ reduction increases in all groups of SBR.

Experimental Validation
However, the experimental and simulated results presented in Figures 4 and 13 display a difference of five orders of magnitude in their G′ values; nevertheless, the overall trends remain consistent. This means that both experimental and simulated outcomes confirm that vulcanized SBR exhibits a distinct Payne effect. The significant difference in magnitude is due to the scale effect of MD simulation. Consequently, the experimental results can provide more reliable evidence to support our simulation approach.

Conclusions
In this study, molecular dynamics (MD) simulations are conducted to study the SBR  Figure 13 shows an increase in the Payne effect in all groups of SBR with an increase in γ 0 . As observed from combining these findings with the results presented in Figure 5, it is apparent that higher S content enhances the deformation resistance of SBR, which is demonstrated by a noticeable upward trend in G ( Figure 13). Furthermore, as γ 0 exceeds 0.1, G reduction increases in all groups of SBR.
However, the experimental and simulated results presented in Figures 4 and 13 display a difference of five orders of magnitude in their G values; nevertheless, the overall trends remain consistent. This means that both experimental and simulated outcomes confirm that vulcanized SBR exhibits a distinct Payne effect. The significant difference in magnitude